module ellipticintegrals_module
! ACM ALGORITHM 577
!
! ALGORITHMS FOR INCOMPLETE ELLIPTIC INTEGRALS
    !
    ! BY B.C. CARLSON AND E.M. NOTIS
    !
    ! ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, SEPTEMBER, 1981.
    !
    !
    !     THIS FILE CONTAINS FOUR SUBROUTINES FOR COMPUTING INCOMPLETE
    !     ELLIPTIC INTEGRALS, FOLLOWED BY SIX DRIVERS FOR TESTING THE
    !     SUBROUTINES.  EACH SUBROUTINE AND EACH DRIVER IS PRECEDED BY
    !     A COMMENT CARD WITH A LINE OF DOLLAR SIGNS, AND EACH DRIVER IS
    !     FOLLOWED BY ITS INPUT DATA IF ANY.  THE FOUR SUBROUTINES HAVE
    !     THE NAMES RC, RF, RD, RJ IN THAT ORDER.  THE DRIVERS HAVE NO
    !     NAMES BUT BEGIN WITH DESCRIPTIVE COMMENTS.  THE FIRST FOUR
    !     DRIVERS TEST RC, RF, RD, RJ IN THAT ORDER.  THE FIFTH DRIVER
    !     TESTS RC AGAINST LIBRARY ROUTINES.  THE SIXTH DRIVER TESTS RF
    !     AGAINST THE FUNPACK SUBROUTINE DELIKM.
    !
    ! NOTE by J. Brandao
    !     This file was lightly reformated, the drivers above mentionned
    !     have been removed and the functions have been encapsulated into a
    !     module.



!interface
!    DOUBLE PRECISION FUNCTION RC(X,Y,ERRTOL,IERR)
!        DOUBLE PRECISION X,Y,ERRTOL
!        INTEGER IERR
!        END FUNCTION RC
!    
!    DOUBLE PRECISION FUNCTION RF(X,Y,Z,ERRTOL,IERR)
!        DOUBLE PRECISION X,Y,Z,ERRTOL
!        INTEGER IERR
!        END FUNCTION RF
!    
!    DOUBLE PRECISION FUNCTION RD(X,Y,Z,ERRTOL,IERR)
!        DOUBLE PRECISION X,Y,Z,ERRTOL
!        INTEGER IERR
!        END FUNCTION RD
!    
!    DOUBLE PRECISION FUNCTION RJ(X,Y,Z,P,ERRTOL,IERR)
!        DOUBLE PRECISION X,Y,Z,P,ERRTOL
!        END FUNCTION RJ
!end interface




contains




DOUBLE PRECISION FUNCTION RC(X,Y,ERRTOL,IERR)
    !
    !          THIS FUNCTION SUBROUTINE COMPUTES THE ELEMENTARY INTEGRAL
    !          RC(X,Y) = INTEGRAL FROM ZERO TO INFINITY OF
    !
    !                              -1/2     -1
    !                    (1/2)(T+X)    (T+Y)  DT,
    !
    !          WHERE X IS NONNEGATIVE AND Y IS POSITIVE.  THE DUPLICATION
    !          THEOREM IS ITERATED UNTIL THE VARIABLES ARE NEARLY EQUAL,
    !          AND THE FUNCTION IS THEN EXPANDED IN TAYLOR SERIES TO FIFTH
    !          ORDER.  LOGARITHMIC, INVERSE CIRCULAR, AND INVERSE HYPER-
    !          BOLIC FUNCTIONS CAN BE EXPRESSED IN TERMS OF RC.  REFERENCE:
    !          B. C. CARLSON, COMPUTING ELLIPTIC INTEGRALS BY DUPLICATION,
    !          NUMER. MATH. 33 (1979), 1-16.  CODED BY B. C. CARLSON AND
    !          ELAINE M. NOTIS, AMES LABORATORY-DOE, IOWA STATE UNIVERSITY,
    !          AMES, IOWA 50011.  MARCH 1, 1980.
    !
    !          CHECK BY ADDITION THEOREM: RC(X,X+Z) + RC(Y,Y+Z) = RC(0,Z),
    !          WHERE X, Y, AND Z ARE POSITIVE AND  X * Y = Z * Z.
    !
    INTEGER IERR,PRINTR
    DOUBLE PRECISION C1,C2,ERRTOL,LAMDA,LOLIM
    DOUBLE PRECISION MU,S,SN,UPLIM,X,XN,Y,YN
    !          INTRINSIC FUNCTIONS USED: DABS,DMAX1,DSQRT
    !
    !          PRINTR IS THE UNIT NUMBER OF THE PRINTER.
    DATA PRINTR/6/
    !
    !          LOLIM AND UPLIM DETERMINE THE RANGE OF VALID ARGUMENTS.
    !          LOLIM IS NOT LESS THAN THE MACHINE MINIMUM MULTIPLIED BY 5.
    !          UPLIM IS NOT GREATER THAN THE MACHINE MAXIMUM DIVIDED BY 5.
    !
    !          ACCEPTABLE VALUES FOR:   LOLIM      UPLIM
    !          IBM 360/370 SERIES   :   3.D-78     1.D+75
    !          CDC 6000/7000 SERIES :   1.D-292    1.D+321
    !          UNIVAC 1100 SERIES   :   1.D-307    1.D+307
    !
    !          WARNING: IF THIS PROGRAM IS CONVERTED TO SINGLE PRECISION,
    !          THE VALUES FOR THE UNIVAC 1100 SERIES SHOULD BE CHANGED TO
    !          LOLIM = 1.E-37 AND UPLIM = 1.E+37 BECAUSE THE MACHINE
    !          EXTREMA CHANGE WITH THE PRECISION.
    !
    DATA LOLIM/3.D-78/, UPLIM/1.D+75/
    !
    !          ON INPUT:
    !
    !          X AND Y ARE THE VARIABLES IN THE INTEGRAL RC(X,Y).
    !
    !          ERRTOL IS SET TO THE DESIRED ERROR TOLERANCE.
    !          RELATIVE ERROR DUE TO TRUNCATION IS LESS THAN
    !          16 * ERRTOL ** 6 / (1 - 2 * ERRTOL).
    !
    !          SAMPLE CHOICES:  ERRTOL   RELATIVE TRUNCATION
    !                                    ERROR LESS THAN
    !                           1.D-3    2.D-17
    !                           3.D-3    2.D-14
    !                           1.D-2    2.D-11
    !                           3.D-2    2.D-8
    !                           1.D-1    2.D-5
    !
    !          ON OUTPUT:
    !
    !          X, Y, AND ERRTOL ARE UNALTERED.
    !
    !          IERR IS THE RETURN ERROR CODE:
    !               IERR = 0 FOR NORMAL COMPLETION OF THE SUBROUTINE,
    !               IERR = 1 FOR ABNORMAL TERMINATION.
    !
    !          ********************************************************
    !          WARNING: CHANGES IN THE PROGRAM MAY IMPROVE SPEED AT THE
    !          EXPENSE OF ROBUSTNESS.
    !
      IF (X .LT. 0.D0  .OR.  Y .LE. 0.D0) GO TO 100
      IF ((X + Y) .LT. LOLIM) GO TO 100
      IF (DMAX1(X,Y) .LE. UPLIM) GO TO 112
  100 WRITE(PRINTR,104)
  104 FORMAT(1H0,42H*** ERROR - INVALID ARGUMENTS PASSED TO RC)
      WRITE(PRINTR,108) X,Y
  108 FORMAT(1H ,4HX = ,D23.16,4X,4HY = ,D23.16)
      IERR = 1
      GO TO 124
  112 IERR = 0
      XN = X
      YN = Y
  116 MU = (XN + YN + YN) / 3.D0
      SN = (YN + MU) / MU - 2.D0
      IF (DABS(SN) .LT. ERRTOL) GO TO 120
      LAMDA = 2.D0 * DSQRT(XN) * DSQRT(YN) + YN
      XN = (XN + LAMDA) * 0.25D0
      YN = (YN + LAMDA) * 0.25D0
      GO TO 116
  120 C1 = 1.D0 / 7.D0
      C2 = 9.D0 / 22.D0
      S = SN * SN * (0.3D0 + SN * (C1 + SN * (0.375D0 + SN * C2)))
      RC = (1.D0 + S) / DSQRT(MU)
  124 RETURN
      END FUNCTION RC




DOUBLE PRECISION FUNCTION RF(X,Y,Z,ERRTOL,IERR)
    !
    !          THIS FUNCTION SUBROUTINE COMPUTES THE INCOMPLETE ELLIPTIC
    !          INTEGRAL OF THE FIRST KIND,
    !          RF(X,Y,Z) = INTEGRAL FROM ZERO TO INFINITY OF
    !
    !                                -1/2     -1/2     -1/2
    !                      (1/2)(T+X)    (T+Y)    (T+Z)    DT,
    !
    !          WHERE X, Y, AND Z ARE NONNEGATIVE AND AT MOST ONE OF THEM
    !          IS ZERO.  IF ONE OF THEM IS ZERO, THE INTEGRAL IS COMPLETE.
    !          THE DUPLICATION THEOREM IS ITERATED UNTIL THE VARIABLES ARE
    !          NEARLY EQUAL, AND THE FUNCTION IS THEN EXPANDED IN TAYLOR
    !          SERIES TO FIFTH ORDER.  REFERENCE: B. C. CARLSON, COMPUTING
    !          ELLIPTIC INTEGRALS BY DUPLICATION, NUMER. MATH. 33 (1979),
    !          1-16.  CODED BY B. C. CARLSON AND ELAINE M. NOTIS, AMES
    !          LABORATORY-DOE, IOWA STATE UNIVERSITY, AMES, IOWA 50011.
    !          MARCH 1, 1980.
    !
    !          CHECK BY ADDITION THEOREM: RF(X,X+Z,X+W) + RF(Y,Y+Z,Y+W)
    !          = RF(0,Z,W), WHERE X,Y,Z,W ARE POSITIVE AND X * Y = Z * W.
    !
    INTEGER IERR,PRINTR
    DOUBLE PRECISION C1,C2,C3,E2,E3,EPSLON,ERRTOL,LAMDA
    DOUBLE PRECISION LOLIM,MU,S,UPLIM,X,XN,XNDEV,XNROOT
    DOUBLE PRECISION Y,YN,YNDEV,YNROOT,Z,ZN,ZNDEV,ZNROOT
    !          INTRINSIC FUNCTIONS USED: DABS,DMAX1,DMIN1,DSQRT
    !
    !          PRINTR IS THE UNIT NUMBER OF THE PRINTER.
    DATA PRINTR/6/
    !
    !          LOLIM AND UPLIM DETERMINE THE RANGE OF VALID ARGUMENTS.
    !          LOLIM IS NOT LESS THAN THE MACHINE MINIMUM MULTIPLIED BY 5.
    !          UPLIM IS NOT GREATER THAN THE MACHINE MAXIMUM DIVIDED BY 5.
    !
    !          ACCEPTABLE VALUES FOR:   LOLIM      UPLIM
    !          IBM 360/370 SERIES   :   3.D-78     1.D+75
    !          CDC 6000/7000 SERIES :   1.D-292    1.D+321
    !          UNIVAC 1100 SERIES   :   1.D-307    1.D+307
    !
    !          WARNING: IF THIS PROGRAM IS CONVERTED TO SINGLE PRECISION,
    !          THE VALUES FOR THE UNIVAC 1100 SERIES SHOULD BE CHANGED TO
    !          LOLIM = 1.E-37 AND UPLIM = 1.E+37 BECAUSE THE MACHINE
    !          EXTREMA CHANGE WITH THE PRECISION.
    !
    DATA LOLIM/3.D-78/, UPLIM/1.D+75/
    !
    !          ON INPUT:
    !
    !          X, Y, AND Z ARE THE VARIABLES IN THE INTEGRAL RF(X,Y,Z).
    !
    !          ERRTOL IS SET TO THE DESIRED ERROR TOLERANCE.
    !          RELATIVE ERROR DUE TO TRUNCATION IS LESS THAN
    !          ERRTOL ** 6 / (4 * (1 - ERRTOL)).
    !
    !          SAMPLE CHOICES:  ERRTOL   RELATIVE TRUNCATION
    !                                    ERROR LESS THAN
    !                           1.D-3    3.D-19
    !                           3.D-3    2.D-16
    !                           1.D-2    3.D-13
    !                           3.D-2    2.D-10
    !                           1.D-1    3.D-7
    !
    !          ON OUTPUT:
    !
    !          X, Y, Z, AND ERRTOL ARE UNALTERED.
    !
    !          IERR IS THE RETURN ERROR CODE:
    !               IERR = 0 FOR NORMAL COMPLETION OF THE SUBROUTINE,
    !               IERR = 1 FOR ABNORMAL TERMINATION.
    !
    !          ********************************************************
    !          WARNING: CHANGES IN THE PROGRAM MAY IMPROVE SPEED AT THE
    !          EXPENSE OF ROBUSTNESS.
    !
      IF (DMIN1(X,Y,Z) .LT. 0.D0) GO TO 100
      IF (DMIN1(X+Y,X+Z,Y+Z) .LT. LOLIM) GO TO 100
      IF (DMAX1(X,Y,Z) .LE. UPLIM) GO TO 112
  100 WRITE(PRINTR,104)
  104 FORMAT(1H0,42H*** ERROR - INVALID ARGUMENTS PASSED TO RF)
      WRITE(PRINTR,108) X,Y,Z
  108 FORMAT(1H ,4HX = ,D23.16,4X,4HY = ,D23.16,4X,4HZ = ,D23.16)
      IERR = 1
      GO TO 124
  112 IERR = 0
      XN = X
      YN = Y
      ZN = Z
  116 MU = (XN + YN + ZN) / 3.D0
      XNDEV = 2.D0 - (MU + XN) / MU
      YNDEV = 2.D0 - (MU + YN) / MU
      ZNDEV = 2.D0 - (MU + ZN) / MU
      EPSLON = DMAX1(DABS(XNDEV),DABS(YNDEV),DABS(ZNDEV))
      IF (EPSLON .LT. ERRTOL) GO TO 120
      XNROOT = DSQRT(XN)
      YNROOT = DSQRT(YN)
      ZNROOT = DSQRT(ZN)
      LAMDA = XNROOT * (YNROOT + ZNROOT) + YNROOT * ZNROOT
      XN = (XN + LAMDA) * 0.25D0
      YN = (YN + LAMDA) * 0.25D0
      ZN = (ZN + LAMDA) * 0.25D0
      GO TO 116
  120 C1 = 1.D0 / 24.D0
      C2 = 3.D0 / 44.D0
      C3 = 1.D0 / 14.D0
      E2 = XNDEV * YNDEV - ZNDEV * ZNDEV
      E3 = XNDEV * YNDEV * ZNDEV
      S = 1.D0 + (C1 * E2 - 0.1D0 - C2 * E3) * E2 + C3 * E3
      RF = S / DSQRT(MU)
  124 RETURN
      END FUNCTION RF




DOUBLE PRECISION FUNCTION RD(X,Y,Z,ERRTOL,IERR)
    !
    !          THIS FUNCTION SUBROUTINE COMPUTES AN INCOMPLETE ELLIPTIC
    !          INTEGRAL OF THE SECOND KIND,
    !          RD(X,Y,Z) = INTEGRAL FROM ZERO TO INFINITY OF
    !
    !                                -1/2     -1/2     -3/2
    !                      (3/2)(T+X)    (T+Y)    (T+Z)    DT,
    !
    !          WHERE X AND Y ARE NONNEGATIVE, X + Y IS POSITIVE, AND Z IS
    !          POSITIVE.  IF X OR Y IS ZERO, THE INTEGRAL IS COMPLETE.
    !          THE DUPLICATION THEOREM IS ITERATED UNTIL THE VARIABLES ARE
    !          NEARLY EQUAL, AND THE FUNCTION IS THEN EXPANDED IN TAYLOR
    !          SERIES TO FIFTH ORDER.  REFERENCE: B. C. CARLSON, COMPUTING
    !          ELLIPTIC INTEGRALS BY DUPLICATION, NUMER. MATH. 33 (1979),
    !          1-16.  CODED BY B. C. CARLSON AND ELAINE M. NOTIS, AMES
    !          LABORATORY-DOE, IOWA STATE UNIVERSITY, AMES, IOWA 50011.
    !          MARCH 1, 1980..
    !
    !          CHECK: RD(X,Y,Z) + RD(Y,Z,X) + RD(Z,X,Y)
    !          = 3 / DSQRT(X * Y * Z), WHERE X, Y, AND Z ARE POSITIVE.
    !
    INTEGER IERR,PRINTR
    DOUBLE PRECISION C1,C2,C3,C4,EA,EB,EC,ED,EF,EPSLON,ERRTOL,LAMDA
    DOUBLE PRECISION LOLIM,MU,POWER4,SIGMA,S1,S2,UPLIM,X,XN,XNDEV
    DOUBLE PRECISION XNROOT,Y,YN,YNDEV,YNROOT,Z,ZN,ZNDEV,ZNROOT
    !          INTRINSIC FUNCTIONS USED: DABS,DMAX1,DMIN1,DSQRT
    !
    !          PRINTR IS THE UNIT NUMBER OF THE PRINTER.
    DATA PRINTR/6/
    !
    !          LOLIM AND UPLIM DETERMINE THE RANGE OF VALID ARGUMENTS.
    !          LOLIM IS NOT LESS THAN 2 / (MACHINE MAXIMUM) ** (2/3).
    !          UPLIM IS NOT GREATER THAN (0.1 * ERRTOL / MACHINE
    !          MINIMUM) ** (2/3), WHERE ERRTOL IS DESCRIBED BELOW.
    !          IN THE FOLLOWING TABLE IT IS ASSUMED THAT ERRTOL WILL
    !          NEVER BE CHOSEN SMALLER THAN 1.D-5.
    !
    !          ACCEPTABLE VALUES FOR:   LOLIM      UPLIM
    !          IBM 360/370 SERIES   :   6.D-51     1.D+48
    !          CDC 6000/7000 SERIES :   5.D-215    2.D+191
    !          UNIVAC 1100 SERIES   :   1.D-205    2.D+201
    !
    !          WARNING: IF THIS PROGRAM IS CONVERTED TO SINGLE PRECISION,
    !          THE VALUES FOR THE UNIVAC 1100 SERIES SHOULD BE CHANGED TO
    !          LOLIM = 1.E-25 AND UPLIM = 2.E+21 BECAUSE THE MACHINE
    !          EXTREMA CHANGE WITH THE PRECISION.
    !
    DATA LOLIM/6.D-51/, UPLIM/1.D+48/
    !
    !          ON INPUT:
    !
    !          X, Y, AND Z ARE THE VARIABLES IN THE INTEGRAL RD(X,Y,Z).
    !
    !          ERRTOL IS SET TO THE DESIRED ERROR TOLERANCE.
    !          RELATIVE ERROR DUE TO TRUNCATION IS LESS THAN
    !          3 * ERRTOL ** 6 / (1-ERRTOL) ** 3/2.
    !
    !          SAMPLE CHOICES:  ERRTOL   RELATIVE TRUNCATION
    !                                    ERROR LESS THAN
    !                           1.D-3    4.D-18
    !                           3.D-3    3.D-15
    !                           1.D-2    4.D-12
    !                           3.D-2    3.D-9
    !                           1.D-1    4.D-6
    !
    !          ON OUTPUT:
    !
    !          X, Y, Z, AND ERRTOL ARE UNALTERED.
    !
    !          IERR IS THE RETURN ERROR CODE:
    !               IERR = 0 FOR NORMAL COMPLETION OF THE SUBROUTINE,
    !               IERR = 1 FOR ABNORMAL TERMINATION.
    !
    !          ********************************************************
    !          WARNING: CHANGES IN THE PROGRAM MAY IMPROVE SPEED AT THE
    !          EXPENSE OF ROBUSTNESS.
    !
      IF (DMIN1(X,Y) .LT. 0.D0) GO TO 100
      IF (DMIN1(X+Y,Z) .LT. LOLIM) GO TO 100
      IF (DMAX1(X,Y,Z) .LE. UPLIM) GO TO 112
  100 WRITE(PRINTR,104)
  104 FORMAT(1H0,42H*** ERROR - INVALID ARGUMENTS PASSED TO RD)
      WRITE(PRINTR,108) X,Y,Z
  108 FORMAT(1H ,4HX = ,D23.16,4X,4HY = ,D23.16,4X,4HZ = ,D23.16)
      IERR = 1
      GO TO 124
  112 IERR = 0
      XN = X
      YN = Y
      ZN = Z
      SIGMA = 0.D0
      POWER4 = 1.D0
  116 MU = (XN + YN + 3.D0 * ZN) * 0.2D0
      XNDEV = (MU - XN) / MU
      YNDEV = (MU - YN) / MU
      ZNDEV = (MU - ZN) / MU
      EPSLON = DMAX1(DABS(XNDEV),DABS(YNDEV),DABS(ZNDEV))
      IF (EPSLON .LT. ERRTOL) GO TO 120
      XNROOT = DSQRT(XN)
      YNROOT = DSQRT(YN)
      ZNROOT = DSQRT(ZN)
      LAMDA = XNROOT * (YNROOT + ZNROOT) + YNROOT * ZNROOT
      SIGMA = SIGMA + POWER4 / (ZNROOT * (ZN + LAMDA))
      POWER4 = POWER4 * 0.25D0
      XN = (XN + LAMDA) * 0.25D0
      YN = (YN + LAMDA) * 0.25D0
      ZN = (ZN + LAMDA) * 0.25D0
      GO TO 116
  120 C1 = 3.D0 / 14.D0
      C2 = 1.D0 / 6.D0
      C3 = 9.D0 / 22.D0
      C4 = 3.D0 / 26.D0
      EA = XNDEV * YNDEV
      EB = ZNDEV * ZNDEV
      EC = EA - EB
      ED = EA - 6.D0 * EB
      EF = ED + EC + EC
      S1 = ED * (- C1 + 0.25D0 * C3 * ED - 1.5D0 * C4 * ZNDEV * EF)
      S2 = ZNDEV * (C2 * EF + ZNDEV * (- C3 * EC + ZNDEV * C4 * EA))
      RD = 3.D0 * SIGMA + POWER4 * (1.D0 + S1 + S2) / (MU * DSQRT(MU))
  124 RETURN
      END FUNCTION RD




DOUBLE PRECISION FUNCTION RJ(X,Y,Z,P,ERRTOL,IERR)
    !
    !          THIS FUNCTION SUBROUTINE COMPUTES AN INCOMPLETE ELLIPTIC
    !          INTEGRAL OF THE THIRD KIND,
    !          RJ(X,Y,Z,P) = INTEGRAL FROM ZERO TO INFINITY OF
    !
    !                                  -1/2     -1/2     -1/2     -1
    !                        (3/2)(T+X)    (T+Y)    (T+Z)    (T+P)  DT,
    !
    !          WHERE X, Y, AND Z ARE NONNEGATIVE, AT MOST ONE OF THEM IS
    !          ZERO, AND P IS POSITIVE.  IF X OR Y OR Z IS ZERO, THE
    !          INTEGRAL IS COMPLETE.  THE DUPLICATION THEOREM IS ITERATED
    !          UNTIL THE VARIABLES ARE NEARLY EQUAL, AND THE FUNCTION IS
    !          THEN EXPANDED IN TAYLOR SERIES TO FIFTH ORDER.  REFERENCE:
    !          B. C. CARLSON, COMPUTING ELLIPTIC INTEGRALS BY DUPLICATION,
    !          NUMER. MATH. 33 (1979), 1-16.  CODED BY B. C. CARLSON AND
    !          ELAINE M. NOTIS, AMES LABORATORY-DOE, IOWA STATE UNIVERSITY,
    !          AMES, IOWA 50011.  MARCH 1, 1980.
    !
    !          CHECK BY ADDITION THEOREM: RJ(X,X+Z,X+W,X+P)
    !          + RJ(Y,Y+Z,Y+W,Y+P) + (A-B) * RJ(A,B,B,A) + 3 / DSQRT(A)
    !          = RJ(0,Z,W,P), WHERE X,Y,Z,W,P ARE POSITIVE AND X * Y
    !          = Z * W,  A = P * P * (X+Y+Z+W),  B = P * (P+X) * (P+Y),
    !          AND B - A = P * (P-Z) * (P-W).  THE SUM OF THE THIRD AND
    !          FOURTH TERMS ON THE LEFT SIDE IS 3 * RC(A,B).
    !
    INTEGER IERR,PRINTR
    DOUBLE PRECISION ALFA,BETA,C1,C2,C3,C4,EA,EB,EC,E2,E3
    DOUBLE PRECISION EPSLON,ERRTOL,ETOLRC,LAMDA,LOLIM,MU,P,PN,PNDEV
    DOUBLE PRECISION POWER4,SIGMA,S1,S2,S3,UPLIM,X,XN,XNDEV
    DOUBLE PRECISION XNROOT,Y,YN,YNDEV,YNROOT,Z,ZN,ZNDEV,ZNROOT
    !          INTRINSIC FUNCTIONS USED: DABS,DMAX1,DMIN1,DSQRT
    !
    !          PRINTR IS THE UNIT NUMBER OF THE PRINTER.
    DATA PRINTR/6/
    !
    !          RC IS A FUNCTION COMPUTED BY AN EXTERNAL SUBROUTINE.
    !
    !          LOLIM AND UPLIM DETERMINE THE RANGE OF VALID ARGUMENTS.
    !          LOLIM IS NOT LESS THAN THE CUBE ROOT OF THE VALUE
    !          OF LOLIM USED IN THE SUBROUTINE FOR RC.
    !          UPLIM IS NOT GREATER THAN 0.3 TIMES THE CUBE ROOT OF
    !          THE VALUE OF UPLIM USED IN THE SUBROUTINE FOR RC.
    !
    !          ACCEPTABLE VALUES FOR:   LOLIM      UPLIM
    !          IBM 360/370 SERIES   :   2.D-26     3.D+24
    !          CDC 6000/7000 SERIES :   5.D-98     3.D+106
    !          UNIVAC 1100 SERIES   :   5.D-103    6.D+101
    !
    !          WARNING: IF THIS PROGRAM IS CONVERTED TO SINGLE PRECISION,
    !          THE VALUES FOR THE UNIVAC 1100 SERIES SHOULD BE CHANGED TO
    !          LOLIM = 5.E-13 AND UPLIM = 6.E+11 BECAUSE THE MACHINE
    !          EXTREMA CHANGE WITH THE PRECISION.
    !
    DATA LOLIM/2.D-26/, UPLIM/3.D+24/
    !
    !          ON INPUT:
    !
    !          X, Y, Z, AND P ARE THE VARIABLES IN THE INTEGRAL RJ(X,Y,Z,P).
    !
    !          ERRTOL IS SET TO THE DESIRED ERROR TOLERANCE.
    !          RELATIVE ERROR DUE TO TRUNCATION OF THE SERIES FOR RJ
    !          IS LESS THAN 3 * ERRTOL ** 6 / (1 - ERRTOL) ** 3/2.
    !          AN ERROR TOLERANCE (ETOLRC) WILL BE PASSED TO THE SUBROUTINE
    !          FOR RC TO MAKE THE TRUNCATION ERROR FOR RC LESS THAN FOR RJ.
    !
    !          SAMPLE CHOICES:  ERRTOL   RELATIVE TRUNCATION
    !                                    ERROR LESS THAN
    !                           1.D-3    4.D-18
    !                           3.D-3    3.D-15
    !                           1.D-2    4.D-12
    !                           3.D-2    3.D-9
    !                           1.D-1    4.D-6
    !
    !          ON OUTPUT:
    !
    !          X, Y, Z, P, AND ERRTOL ARE UNALTERED.
    !
    !          IERR IS THE RETURN ERROR CODE:
    !               IERR = 0 FOR NORMAL COMPLETION OF THE SUBROUTINE,
    !               IERR = 1 FOR ABNORMAL TERMINATION.
    !
    !          ********************************************************
    !          WARNING: CHANGES IN THE PROGRAM MAY IMPROVE SPEED AT THE
    !          EXPENSE OF ROBUSTNESS.
    !
      IF (DMIN1(X,Y,Z) .LT. 0.D0) GO TO 100
      IF (DMIN1(X+Y,X+Z,Y+Z,P) .LT. LOLIM) GO TO 100
      IF (DMAX1(X,Y,Z,P) .LE. UPLIM) GO TO 112
  100 WRITE(PRINTR,104)
  104 FORMAT(1H0,42H*** ERROR - INVALID ARGUMENTS PASSED TO RJ)
      WRITE(PRINTR,108) X,Y,Z,P
  108 FORMAT(1H ,4HX = ,D23.16,4X,4HY = ,D23.16,4X,4HZ = ,D23.16, 4X,4HP = ,D23.16)
      IERR = 1
      GO TO 124
  112 IERR = 0
      XN = X
      YN = Y
      ZN = Z
      PN = P
      SIGMA = 0.D0
      POWER4 = 1.D0
      ETOLRC = 0.5D0 * ERRTOL
  116 MU = (XN + YN + ZN + PN + PN) * 0.2D0
      XNDEV = (MU - XN) / MU
      YNDEV = (MU - YN) / MU
      ZNDEV = (MU - ZN) / MU
      PNDEV = (MU - PN) / MU
      EPSLON = DMAX1(DABS(XNDEV),DABS(YNDEV),DABS(ZNDEV),DABS(PNDEV))
      IF (EPSLON .LT. ERRTOL) GO TO 120
      XNROOT = DSQRT(XN)
      YNROOT = DSQRT(YN)
      ZNROOT = DSQRT(ZN)
      LAMDA = XNROOT * (YNROOT + ZNROOT) + YNROOT * ZNROOT
      ALFA = PN * (XNROOT + YNROOT + ZNROOT) + XNROOT * YNROOT * ZNROOT
      ALFA = ALFA * ALFA
      BETA = PN * (PN + LAMDA) * (PN + LAMDA)
      SIGMA = SIGMA + POWER4 * RC(ALFA,BETA,ETOLRC,IERR)
      IF (IERR .NE. 0) GO TO 124
      POWER4 = POWER4 * 0.25D0
      XN = (XN + LAMDA) * 0.25D0
      YN = (YN + LAMDA) * 0.25D0
      ZN = (ZN + LAMDA) * 0.25D0
      PN = (PN + LAMDA) * 0.25D0
      GO TO 116
  120 C1 = 3.D0 / 14.D0
      C2 = 1.D0 / 3.D0
      C3 = 3.D0 / 22.D0
      C4 = 3.D0 / 26.D0
      EA = XNDEV * (YNDEV + ZNDEV) + YNDEV * ZNDEV
      EB = XNDEV * YNDEV * ZNDEV
      EC = PNDEV * PNDEV
      E2 = EA - 3.D0 * EC
      E3 = EB + 2.D0 * PNDEV * (EA - EC)
      S1 = 1.D0 + E2 * (- C1 + 0.75D0 * C3 * E2 - 1.5D0 * C4 * E3)
      S2 = EB * (0.5D0 * C2 + PNDEV * (- C3 - C3 + PNDEV * C4))
      S3 = PNDEV * EA * (C2 - PNDEV * C3) - C2 * PNDEV * EC
      RJ = 3.D0 * SIGMA + POWER4 * (S1 + S2 + S3) / (MU * DSQRT(MU))
  124 RETURN
      END FUNCTION RJ

end module ellipticintegrals_module
